%% Scripts to read all the data and select a sub set in time and save the data to a file
% wrote to extract data for storm time induction
% Zurich July 11, 2013
%close all;
clear

% Get the directory for netcdf database
%localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/absolute/';
localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/secular_removed/';
%localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/Canada_Absolute_netCDF/';
%localdirectory = '/nfs/satmag_work/mnair/projects/obs_mag_data/update_2012/Image_Absolute_netCDF/';
%localdirectory = '/Users/manojnair/data/obs_mag_data/Canada/Absolute_netCDF/';
%localdirectory = '/Users/manojnair/data/obs_mag_data/image/Absolute_netCDF/';
%localdirectory = '/Users/manojnair/data/obs_mag_data/Greenland/Secular_Removed/';
% Get a listing of the files
S = dir([localdirectory '*.nc']);
for i= 1:length(S),
        filestring(i,:) = S(i).name;
end;

% Time array for the data
% DO NOT CHANGE THIS unless the netCDF file changes its length
time_array_min = (datenum(1995,01,01,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a= length(time_array_min);
where_to_store = 0;
length_of_record = a;
latitude_max = 55.0; % latitude cut off (abs)



% start date for which data is required
%start_date = datenum(2011,10,20);
%start_date = datenum(2003,10,26, 12, 0, 0);
start_date = datenum(2003,11,09, 0, 0, 0);
%start_date = datenum(2009,07,17);
%start_date = datenum(2011,08,1);
%start_date = datenum(2012,07,10);
start_index = round((start_date - time_array_min(1)) * 1440) + 1;


% length of the extracted data 
n_length = 10 * 1440;
xaxis = 1:n_length;

%number of available stations
n_available = 1;

for i = 1:length(S),

nc_fname = [localdirectory S(i).name];

fprintf('%d %s\n',i,nc_fname);   

S1 = dir(nc_fname);

if ~isempty(S1),

% read the data for a station    
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] =  ... 
    read_geomag_netcdf(nc_fname, start_index, n_length, 1);
end;

% select stations only for given latitude bands. remove stations with large
% data gaps

% In some USGS observatories, the lat/long are given as characters
if ischar(obj.geospatial_lat)
    obj.geospatial_lat = str2double(obj.geospatial_lat);
    obj.geospatial_lon = str2double(obj.geospatial_lon);
end;


% Get the geomagnetic latitude
tg = 90 - obj.geospatial_lat;
pg = obj.geospatial_lon;

[tm, pm] = gg2gm_new(tg,	pg);
maglat = - (tm - 90);


if abs(maglat) <latitude_max & sum(isnan(x_data)) < 1440 & ...
        sum(isnan(y_data)) < 1440 & sum(isnan(z_data)) < 1440 ,
    L = isnan(x_data);
    x_data = interp1( xaxis(~L) , x_data(~L), xaxis, 'linear',nanmedian(x_data));
    L = isnan(y_data);
    y_data = interp1( xaxis(~L) , y_data(~L), xaxis, 'linear',nanmedian(y_data));
    L = isnan(z_data);
    z_data = interp1( xaxis(~L) , z_data(~L), xaxis, 'linear',nanmedian(z_data));
    
    if sum(isnan(x_data)) > 0,
        break;
    end;
    
    X(n_available,:) = x_data;
    Y(n_available,:) = y_data;
    Z(n_available,:) = z_data;
    
    name(n_available,:) = S(i).name(1:3);
    theta_gg(n_available) = 90 - obj.geospatial_lat; % co latitude
    phi_gg(n_available) = obj.geospatial_lon;
    theta_gm(n_available) = tm;
    phi_gm(n_available) = pm;
    n_available = n_available + 1;
end;

end;

fday = time_array_min(start_index : start_index + n_length - 1);

%save /nfs/satmag_work/mnair/projects/obs_mag_data/obs_Apr_2000 X Y Z name phi_gg theta_gg fday
%save /Users/manojnair/data/obs_mag_data/obs_Apr_2000_image X Y Z name phi_gg theta_gg fday
%save /Users/manojnair/data/obs_mag_data/obs_Apr_2000_canada X Y Z name phi_gg theta_gg fday

save /nfs/satmag_work/mnair/projects/substorm/obs_NovSq_2003 X Y Z name phi_gg theta_gg theta_gm phi_gm fday

